module general
use omplib
use technical
implicit none
save
!
private
public :: linear_interp_to_offgrid,get_grid_loc,foldback_pbc,peskin_delta,cross,curlF
! ********************************************************* !
contains
!!-----------------------------------------------------------
subroutine linear_interp_to_offgrid(xp,yp,n1,n2,deltax,deltay,psi,psi_interpolated)
!
! Linearly interpolates a quantity to an off-grid point
!
  double precision :: xp,yp
  double precision, intent(in) :: deltax,deltay
  integer :: xgrid,ygrid,xgrid_next,ygrid_next 
  double precision :: ddx,ddy,dxdy,g1,g2,g3,g4
  integer :: n1,n2
  double precision,dimension(n1,n2) :: psi
  double precision :: psi_interpolated
!
  call get_grid_loc(xp,n1,deltax,xgrid,xgrid_next,ddx)
  call get_grid_loc(yp,n2,deltay,ygrid,ygrid_next,ddy)
!
  g1=psi(xgrid,ygrid)
  g2=psi(xgrid_next,ygrid)
  g3=psi(xgrid,ygrid_next)
  g4=psi(xgrid_next,ygrid_next)
!
  psi_interpolated = g1+(ddx/deltax)*(g2-g1)+&
    (ddy/deltay)*(g3-g1) + (ddx*ddy/(deltax*deltay))*(g1+g4-g2-g3)
!
endsubroutine linear_interp_to_offgrid
!!-----------------------------------------------------
!!
subroutine get_grid_loc(xip,Nbox,dxi,xi_grid,xi_gridn,ddxi)
  implicit none
  double precision, intent(in) :: xip,dxi
  integer,intent(in) :: Nbox
  integer,intent(out) :: xi_grid,xi_gridn
  integer :: xiint
  double precision,intent(out) :: ddxi
!--------------------
  xiint=int(xip/dxi)
  ddxi=xip-xiint*dxi
  xi_grid=mod(xiint,Nbox)
  if (xip .lt. 0.) then
    ddxi=dxi+ddxi
    xi_grid=Nbox-1+xi_grid
  endif
  xi_gridn=mod(xi_grid+1,Nbox)
!
! because it is fortran we add one to indices
!
  xi_grid=xi_grid+1
  xi_gridn=xi_gridn+1
end subroutine get_grid_loc
!! --------------------------------------------------
subroutine foldback_pbc(ii,ix,n1)
  integer, intent(out) :: ii
  integer, intent(in) :: ix,n1
  if(ix.gt.0)then
    ii = 1 + mod(ix-1,n1)
  else
    ii = n1 + mod(ix,n1)
  endif
endsubroutine foldback_pbc
!! ---------------------------------------------------
function peskin_delta(rx,h)
  double precision,intent(in) :: rx,h
  double precision :: rr
  double precision :: peskin_delta
  rr = rx/h
  if (rr .gt. 2.) then
    peskin_delta = 0.
  else
    peskin_delta =   (1./h)*(1./4.)*( 1+cos(pi*rr/2.) )
  endif
endfunction peskin_delta
!! ---------------------------------------------------
function cross(A,B)
  double precision, dimension(3),intent(in) :: A,B
  double precision,dimension(3) :: cross
  cross(1) = A(2)*B(3)-A(3)*B(2)
  cross(2) = A(3)*B(1)-A(1)*B(3)
  cross(3) = A(1)*B(2) - A(2)*B(1)
endfunction cross
!! ---------------------------------------------------
function curlF(kk,A)
  double complex, dimension(3),intent(in) :: A
  double precision, dimension(3), intent(in) :: kk
  double complex,dimension(3) :: curlF
  curlF=zi*dcmplx( cross( kk,real(A) ), cross( kk,aimag(A) ) )
endfunction curlF
!! ---------------------------------------------------


! ********************************************************* !
endmodule general
